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CN , Abstract 

l'^ ' The simulation of the time dependent evolution of the resonant tunneling diode is done by a 

multiscale algorithm exploiting the existence of resonant states. After revisiting and improving 
the algorithm developed in [N. Ben Abdallah, O. Pinaud, J. Comp. Phys. 213 (2006) 288- 
310] for the stationary case, the time dependent problem is dealt with. The existence of two 
resonances corresponding to the initial potential and to the local time potential lead to the 
decomposition of the wave function into a non resonant part and two resonant ones. The 

Ph ' resonant parts are dealt with by a projection method. The simulation times are shown to be 

■^T , reduced by a factor two. 
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1 Introduction 

>\ 

^, . Wc arc interested in the simulation of open resonant systems, typically the resonant tunneling 

diode. Because of the small size of these devices, quantum effects are observed and the Schrodinger 

(^ . equation is required for the modeling. Classical observables like electronic and current densities 

are given by integrals on the energy involving generalized eigenfunctions. The resonant tunneling 
diode presents resonances and the variation of generalized eigenfunctions tends to be singular in 

f^*) ' the vicinity of resonant energies. Therefore, to compute correctly the densities, the integral on the 

energy variable needs a refined mesh next to these resonances and a big number of Schrodinger 
equations needs to be solved. Adaptive refinement was developed in |16j in the stationary regime 
for a resonant tunneling diode. The method does however not extend to the transient case since 

/\ ' the resonant energies do move as time varies. Lately, inspired by an idea for the transient case 

from [TTj, an algorithm was proposed by [5] for the discretization of the stationary case. The 
method consists in decomposing the wave function in a part living in the well between the double 
barrier and a second part which is mostly localized outside this well. The latter is called the non 
resonant part and the former is the resonant one. The non resonant part is smooth with respect 
to the energy variable and only requires a coarse energy mesh. At variance, the resonant part has 
sharp peaks at resonant energies. It is computed by a projection method after a precomputation 
of resonant states. In the present work, we first present an improvement of 0, the resonance 
being computed more precisely by solving the non-linear eigenvalue problem verified by the res- 
onance and written in [T3]. Then we present a transient algorithm inspired by the idea of |17) . 
More precisely, the algorithm consists in the computation of the resonance at each time step, the 
resonant information is therefore captured and the frequency mesh can be chozen coarse enough 
thus reducing the simulation time. The resonant part of the wave fonction has two components, 
one which corresponds to the dissipation of the initial resonance and one which is proportional 
to the resonant mode. Let us mention that in [S], a new finite element method using the WKB 
approximation was developed to reduce the numerical cost (see also [21 [13] )• Since we want to 
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concentrate on the energy discretization we will use a standard finite difference method for the 
space discretization. 



2 The model 

The considered device is assumed to occupy the domain [0,L], L > 0. The model consists in 
an infinite number of Schrodinger equations coupled to the Poisson equation. The Schrodinger 
equation involves the time dependent Hamiltonian 



H{t) = -—dl + Uit) + V{t), 
Im 

where h is the reduced Planck constant, m is the effective mass of the electron and x ^ M \s the 

position variable. 

The external potentiel U{t) = Vq + B{t) is given as if a data of the problem where 



the well W is given by 
and the applied bias by 



Vq = I'ol [02,62] + ^' 



oi — ai 



with uq ^ and B{t) > scalars representing respectively the height of the barrier and the 
amplitude of the applied bias in eV^ and we have: 

< ai < a2 < as < 63 < 62 < 61 < L. 

The points ai and 61 are the extremities of the diode. 
Such a potential U{t) is represented in Figure [TJ 
The nonlinear potential V{t) is due to Coulomb interaction and satisfies the Poisson equation: 
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Figure 1: External potential U{t). 



-dlV=^-{n[V]-nn). (0,i) 
y(0) = V{L) = 0, 



(2.1) 



where q is the elementary charge of the electron, e is the dielectric constant, no is the doping equal 
to 

no = ri£,(l[o,ai[ + l]6i,L]) + i^oho-iM]' 

with n\) > njj > 0, and the electron density n[V] is given by: 



n[V]ix)^ g{k)\<^kix)\'dk, xe[0,L], 



(2.2) 



in the stationary regime and 

n[V]{t,x)^ f gik)\^kit,x)\^dk, xe[0,L], (2.3) 



in the time dependent regime. Here g is the injection profile and the $fc, resp. the vE'/c(i), are the 
wave functions corresponding to the Hamiltonian H{0), resp. H{t), for a wave vector k G JR. 
We will choose the function g to be the one dimensional Fermi-Dirac integral : 



2^2ft2 "'^- ' -^-y k^T J J' 

where ks is the Boltzmann constant, T is the temperature of the semiconductor and Ep is the 
Fermi level. 

In the stationary case, for a given initial bias Bj — B{0), wc will note Uj = U{0) the corresponding 
initial exterior potential and Vi a corresponding solution to (|2.ip(j2.2p . then, the wave functions 
are the restriction to [0, L] of the solutions to the stationary Schrodinger equation 

-^^'^k + {Ui + Vi)<^k=Ek<^k, xelR, (2.4) 

Zm ax-' 

with scattering conditions 

f $fc(x) =e*'=^+r(/c)e-*'^^, a; < 

for fc > 



[ $fe(x) = i(fc)eV'='+2™-B,/;i2^^ x> L, 

and 

f ^k{x) = t(fc)e-V'='-2™s,/ft2^^ ^^^ 

\ $fc(x) = e'*^^ + r(fc)e-*'==^, x > L, 

where 



for fc < , 



, „ , fc > 

^^= .IT. _ - ^ (2-5) 



ri^fc^ 



. 2™ Bj, fc<0. 

In the transient case, the wave functions are the restriction to [0, L] of the solutions to the time- 
dependent Schrodinger equation 

i ihdt^k{t) = ~£^dl'$k{t) + {u{t) + v{t))^k{t), xeM 

\ *fe(0) = $fe, 
where V{0) = Vj. 

3 The stationary algorithm revisited 

3.1 Recalling the standard algorithm 

To solve the stationary non- linear Schrodinger- Poisson system p.ip (j2.2|) p.4|) , we will use the same 
algorithm as in J16j . It is based on a Gunimel iteration, see |11] , which consists in the computation 
of a sequence Vj, / > where the potential Vj^^ at step / -f 1 is deduced form the potential V} at 
step / by solving the following non-linear equation: 

f -£.V}+' = ^{n[Vl]expiiV^-VJ+')/Vref)-nn), (0,L) 
1 V}+\0)^VJ+\L)^0, 



for a fixed reference potential Vref- 

The repartition function g being exponentially decreasing at infinity, the integral in (|2.2p can be 
restricted, in computations, to a domain [— A:j\f , fcj\/] where kM is chosen to be big enough. We 
consider a uniform discretization xq = 0, xi, ..., Xj, ..., xj = L of the interval (0, L) with fixed mesh 
size Ax and a (non-necessarly uniform) discretization kg = —kM,ki,...,kp,...,kp ~ ku of the 
interval [— /ca/,^^/]- Then, for a given initial potential Vj", the algorithm writes: 

Algorithm 3.1 (Gummel algorithm). 

Fbc ; = 0. 

Do While IIT//+1 - Y\\\2 > e: 

Computation of the density 

5*1. For p = 0,...,P: computation of $' wave function corresponding to the potential Vj 
and the frequency kp. 

S2. Numerical integration: for j = 0, ..., J 

dk. 



p=0 


where 


4. = 





CoupHng to the Poisson equation : Gummel iteration 

Computation of the potential V}^^ from the potential Vj and the density ti' by solving 
equation (|3.ip . 

Set ; = / + i. 

End Do . 

The equation (|3.ip is nonlinear and it is solved using a Newton method where the Laplacian is 
discretizcd with finite differences. 

Because of the peaked form of the transmition near resonances, the method used to compute the 
density n[VJ] from the potential Vj (steps SI and 5*2) is of major importance, and the Gummel 
algorithm may fail to converge if a non accurate method is used. 

In |16| . the steps SI and 52 are performed without particular treatment of resonances: the wave 
functions are computed on [0, L] solving the Schrodinger equation with transparent boundary 
conditions and the integrals /' are computed with a trapezoidal rule. In that case, the convergence 
is provided by the choice of the frequency mesh {kp, p = 0, ..., P}: the mesh size fcp+i — kp is refined 
when a resonance peak is detected. As said in the introduction, this forces to solve a big number 
of Schrodinger equations and increases the numerical cost. However, this process, that we will 
call Direct Resolution, gives good results and will be used to evaluate the performance of the One 
Mode Approximation algorithm presented in what follows. 

As proposed in [5], the One Mode Approximation consists in decomposing the wave function in 
a non-resonant part and a resonant part proportional to the first resonant mode. In [5], using a 
WKB interpolation to compute each part of the wave function, an adapted treatment of the step 
51 is realized and the convergence of the Gummel method is possible. 

In the present work, we do not use the WKB interpolation, however the accuracy required at step 
SI is reached by an improvement of the computation of the first resonance. In particular, a precise 
computation of its imaginary part is essential. A simple reconstition of the wave functions not 
allowing to modify the frequency mesh, we obtain convergence with a large frequency mesh by 
adapting the step 52 such that all the resonant information is taken into account. It is done by 
an explicit integration of the coefficient of proportionality with the resonant mode instead of the 
trapezoidal rule. 



3.2 Accurate computation of resonances 

In the founder work [T] , see also [H] , resonances of a sel-adjoint operator are defined using analytic 
transformations, and it is a common fact that it correspond to an eigenvalue in a modified L^ space 
[5] . Using the second approach, it is shown in [5] and [13] that the resonances of the Hamiltonian: 

-|;^'+* 

where the potential Q verifies: 

g(x) == 0, x < 0, and Q{x)^Ql,x>L, 
are the complex values z such that 3u G H^{0,L) verifying ||w||£,2/q^) = 1 and: 

-^u'{0) + ts{z)u{0)^0, (3.2) 

Here s(z) denotes the determination of the square root which is holomorphic on C \ ilR- and 
defined as follows: for z = pe'^, with p > and £ (— f , ^), s{z) = ype*2 . This determination is 
different from the one used in [5] and leads to a convergent Newton algorithm as will be illustrated 
herebelow. In [S] , the computation of the resonance is done by an iterative method which converges 
slowly. The homogeneous transparent boundary conditions in the problem above give a restriction 
of the eigenproblcm for the resonant mode to the domain (0, L) and allows to make numerical 
computations. In particular, we introduce a variational formulation for problem p.2p including 
the boundary conditions such that the discretization leads to a non-linear eigenvalue problem 
solved with a Newton-like algorithm. 

Multiplying the first equation in (|3.2|) by the conjugate of a function v G H^{0,L) and integrating 
by part, we get the following equation on the resonant mode: 

h"^ [^ — [^ h 

/ u'v'dx+ Quvdx-i^^[s{z){uv){0) + s{z + QL){uv){L)] 

2m Jo Jo v2to 

z I uvdx = 0, 







where the boundary conditions were considered. The discretization is performed using a finite 
element method with P^ function basis: pj, j = 0,..., J defined by Pj{xj') ~ 5jji. This leads to 
the following problem: find (u, z) G C"^+^ x C such that: 

(3.3) 




where M{z) is a non-linear matrix valued function of z defined by: 

M{z) = Ml + s{z)M2 + s{z + Ql)M3 - zA'U, 

and the matrices Mi,...,M4 G A^j+i(C) are given in Appendix |21 This kind of problem was 
studied in [T3] and [TU]. The term u^u appearing in problem p.3p is not differentiable with 
respect to u, therefore a Newton method must be modified to be used here. Following [15], for a 
given iterate (u",z") verifying (u")^u" = I, we are looking for a direction (5u",(5z") such that 
(u" -|- (5u",z" + (5z") is solution to problem p.3[) . Using (u")^u" = I the system: 

J 7\/(z" + 5z")(u" + 5u") = 
1 (u" + d^u")^(u" -f (5u") = I 



gives at order 2: 



M(z")5u" + (5z"7\/'(z")u" = -M(z")u" 
(u")^(5u" + ((5u")^u" = 0. 



(3.4) 



As remarked earlier, there is a problem with the second equation, however, it is enough to impose 
(u")^5u" = in order to verify p.4p . Therefore, we obtain the following linear system: 



M(z") M'(z")u" 
(u")^ 



(5z" 



— r 




where r" == A/(z")u". 



(3.5) 



Its resolution corresponds to an iteration of our Newton-like method to compute the resonant 
mode. The assumption (u")^u" = 1 that we made to obtain this system is verified at order 2 as 
long as (u")^u" = 1. 



3.3 The One Mode Approximation 

We start with the descritpion of the step 5*1, in the One Mode Approximation, of the computation 
of a wave function $fc solution to the stationary Shrodinger equation (j2.4p for a given frequency 
k. Following the works |5] and |17j . the One Mode Approximation consists in the decomposition: 



<^k = ^l 



*L 



where the non-resonant part $"' solves the stationary Schrodinger equation: 



r ^' d^ rr 



nr 
k 



Ek<i>T, 



X e Ft, 



(3.6) 



with filled potential Uijm = Ui + vol[a3,b3], and where we omit the indice I appearing in Algorithm 
13.11 to simplify notations. The function $^'" is computed on [0,i] by the resolution of p.6p with 
the exact transparent boundary conditions: 

f (*r)'(0)+*fc*r(0) = 2ifc. 



($^'')'(L) - i^k^ + 2mBi/h?^Y{L) = 0, 



for fc > 0, and: 



j {*r)'(0) + W^^ - 2mBi/h?^l'-{Q) = 0, 
\ {<^1'')'{L) + ik<^Y{L) = 2ikeJ''^, 

for k < 0. This resolution is perfomed with a RK4 method as it is done to compute ^k in the 
Direct Resolution, see [16] . The statistic g being fastly decreasing at infinity, the resonant part $^ 
is searched, on (0,L), proportional to the resonant mode uj of mimimal resonant energy Re(z/) 
where [—-^-^ +Ui + V7]m/ = zjuj and L \uj{x)\'^dx = 1. In other words, we look for $^. of the 
form: 

'Pl{x)=9kui{x), xe{0,L). 

Like in [S], the condition that $fc verifies the stationary Schrodinger equation (|2.4p gives the 
following explicit value of the proportionality coefficient: 



1 



(3.7) 



7k = — I'o / ^l^'uidx. 

Zl - -fcfe J as 

The resonance and the resonant mode are computed using the method presented in section [3T2] with 
the potential Q — Uj -\-Vi. The method is initialized at the fondamental energy and fondamental 
mode of the Hamiltonian 



equipped with homogeneous Dirichlet boundary conditions at a2 and 62- It's shown in [5] [7] that 
the real part of resonances are weU approached by the eigenvalues of the Dirichlet Hamiltonian 
p.8[) . The imaginary part of resonances being small, such an initialization garanties that the 
algorithm converges to the resonance with smaller energy. This achieves the step SI. 
For the step S2, the peaked form p.7p of the coefficient of proportionality with the resonant mode 
is used in the computation of the intregral: 



ijp+i 



9ik) i^k), 



dk. 



Using an argument of localisation of support, see |17| . the cross term in the development of 



in 



hui) 



can be neglected, which gives the approximation: 



^pj 



9ik) (*r). 



dk 



gik)\0kfdk\iui),f 



The non resonant part being regular with respect to fc, the first integral can be computed with a 
simple trapezoidal rule. The keypoint is then the approximation of Jp := /t''^^ g{k) \0k\ dk. Using 
p.7p . we have: 

1 



Jr, 



Rk 



\Ek~zi\ 



rdk , 



where 



Rk = g{k)vl 



&3 



^l^'ujdx 



The wave function $5!^ has no resonance peak and therefore, Rk is regular with respect to k. Then, 
we can make the linear interpolation on [fcp, fcp+i]: 



where 



This gives 



R 



■P+i 



Afc 



Rk = ctpk + f3p , 



— ^"^ ^^ = A^ 



Jp — OLp 



— -^dk + Pp / — ^ 

\Ek-zi\ Jkp \Ek~-zi\ 



dk. 



(3.9) 
(3.10) 



Using (|2.5p and writing zj = Ej — i^ with Ej, F/ > 0, a direct computation gives: 



dk = j X (kp,kp+i,ci,di) , 



and 



/fcp \Ek~zi\ 
where the function x" is given in Appendix [C] and 



\Ek-zi\ 



dk = "f X {kp,kp+i,ci,di) 
(3.11) 



2m \ , T/ 

This approximation has the advantage to be pcackcd around the resonance. 

Remark 3.2. In numerical applications, the computation of the resonance is done only one time 
at each step of the integer I, the resonance is the same for all integer p. We note also that the 
integral Jp can be computued only one time at each step of the integer p, it is the same for all 
integer j. 



4 The time-dependent algorithm 

4.1 The algorithm 

As in the stationary regime, we first present the algorithm proposed in |16| to solve the time- 
dependent non-linear Schrodinger-Poisson system (j2.ip(|2.3l)(|2.6p . Then, we give the details of a 
time-dependent one-mode approximation algorithm. 

The space and frequency meshes are defined as in section 13.11 For a given exterior potential 
U{t) and time step At, the algorithm corresponds to the computation of the sequence F' of 
approximations of the self-consistant potential at time i/ = Z At. The initial potential V'^ and wave 
functions '^'2, are given by Vj, and respectively $p, discret solution to the stationary non- linear 
Schrodinger-Poisson system (|2.ip(|2.2p(|2.4p with exterior potential equal to Uj. The resolution of 
the Schrodinger equation ()2.6p involves a Crank- Nicolson scheme which is semi-implicit and insures 
stability. Therefore, the value of the intermediary potential l/'+2 must be computed to perform the 
algorithm. It is done by an extrapolation using the Gummel iteration (|3.ip . Then, the algorithm 
writes: 

Algorithm 4.1 (Transient algorithm). 
Do For / > 0: 

Computation of the intermediary potential 

Computation of the potential T/'+2 from the potential V'and the density n' by the resolution 
of the non-linear equation: 

f -^^'+^ = 4("[^']exp((F'-y'+^)/Ke/) -«!.), (0,L) 

for a fixed reference potential V^e/- 

Computation of the density 

5*1. For p = 0, ..., P: computation of '^i^^ wave function at time t'+^ and frequency kp using 
the potential X^'+2 and the wave function ^f' 



S2. Numerical integration: for j ~ 0, ..., J 

^/;+\ where /^^ = / g{k) 



-r 



'A^'] 



2 
dk. 



Coupling, Poisson equation 

Computation of the potential y'+^ from the density n'+^ by solving the Poisson equation 
(EH): 



.£,ym^^(„/+i_„^)^ (0,L) 



(4.2) 
F'+i(0) = V^'+i(i)=0. 



End Do. 



As in the stationary case, the equation (|4.ip is solved using a Newton method. The linear 
equation (j4.2p is discretized with finite differences and solved by a simple matrix inversion. 
Due to the presence of resonances, the steps 5*1 and 5*2 to compute density will be crucial here 
also. 

In the algorithm proposed in [T^, which we will call equally Direct Resolution, the steps 51 and 
S2 are performed without particular treatment of resonances: the wave functions are computed on 
[0,L] solving the Schrodinger equation using a Crank-Nicolson Scheme with discrete transparent 



boundary conditions and the integrals /' ■ are computed with a trapezoidal rule. In that case, the 
accuracy of the method is provided by the imposition of a uniform frequency mesh with a small 
mesh size everywhere (an other strategy would be to choose a small frequency mesh size only in 
the region where will live the resonant energy). Indeed, if the mesh size is small only near the 
initial resonant energy then, because of the time evolution of the resonance, the refined mesh will 
loose its relevance. Therefore, for the Direct Resolution, the number of Schrodingcr equations to 
solve in the time dependent regime is more important than in the stationary regime. Moreover, 
the Crank-Nicolson method to solve one Schrodingcr equation requires a matrix inversion of big 
size (equal to J + 1) and is numerically much more expensive then the RK4 stationary method. 
In this context, it is important to look for an adapted treatment of the resonance peaks to reduce 
the number of frequency points. 

In the following section, we propose a one-mode approximation method which extends the method 
proposed in section [3.31 to the time-dependent case. 



4.2 The time dependent One Mode Approximation 

Recall first that in the initial work ^7\, the One Mode Approximation was presented in the time 
dependent case with a simplified model. Then, let us start with the descritpion, in the One 
Mode Approximation, of the decomposition of a wave function *fe(t) solution to the transient 
Shrodingcr equation (j2.6p for a given frequency k. The One Mode Approximation consists in the 
decomposition: 

^k{t)^^Tit) + mt). (4.3) 

where the non-resonant part ^^^{t) solves the transient Schrodingcr equation: 



wm = * 



(4.4) 



k 



with filled potential Ufui{t) = U{t) + uolfaa.b.,], where $)!'' is solution to p.6p with the initial 
potential Uijm + Vj. In comparison with the stationary algorithm, it seems natural to look for 
^^(f) proportional to the resonant mode u{t) corresponding to the first resonance z{t) of the 
Hamiltonian H(t) at time t which verifies 

i-^d'. + U{t) + V{t)]u{t) = z{t)u{t) (4.5) 

and /p \u{t,x)\'^dx = 1. However, since the change of potential at time t = is abrupt, the 
adiabaticity hypothesis of [T7] is not satisfied and one expects two peaks at both the energy of the 
resonant mode at time t = 0~ and the resonant energy at time t. This is what happens numerically. 
Indeed, if denoting z{t) = E(t) — i—^, where E{t), T{t) > 0, in Figure[2]are plotted the logarithm 
C(i, k) of the charge in the well for one wave function with respect to the frequency k: 

C(t,fc) -log I J ' \^kit,x)\''dx\ (4.6) 

and the frequencies corresponding to the resonant energy 



fcfl W = - V 7^(^W + Bit)) , k+{t) = \I^E{t) , (4.7) 

at different time t. The resolution is performed with the external potential U{t) corresponding 
to the bias B{t) defined as follows: for the biases Bj = OeV and Boo = 0.1 eF, and for the time 
to = 10-12 g^ ^Q j^a,ve B{t) = Bj for t < 0, B{t) ^ Boo for t > to and on (0,io), B{t) is the 
polynomial of degree 3 such that B{t) is a C^ function on M. The method used is the Direct 
Resolution with the numerical parameters J — 300, P — 1500, Ai = IQ-^^s and the physical 



t=D.6x10 s 





C(t,k) 





t=2x10"^^s 
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t' 
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■ — ' 


C(t,k) 

___k„(l) 
- - k*(l) 






Figure 2: Logarithm C{t, k) of the charge in the well for one wave function with respect to the 
frequency k and frequencies, k~^{t) and k~j^[t), corresponding to the resonant energy, at different 
time t, for a time dependent bias. 
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parameters given in section [5] 

We see that for small time the resonance peak is not localised at the resonant energy E{t) at time 
t but it's localised at the initial resonant energy Ej. Then, the initial peak vanishes while a peak 
grows at E{t). Therefore, we look for ^^ of the form: 

%{t,x)^ekv{t,x) + Xk{t)u{t,x), xe(0,L), (4.8) 

where v(t) is the propagation of the initial resonant mode by the Schrodinger equation: 

( ^h^Mt)^[-i:^l + uit) + v{t)]v{t), xeR, ^^g^ 

1 v{0) = ui 

and Ok is the proportionality coefficient in the stationary regime corresponding the initial data 
such that: 

$fc(x) = $r(^) + Okuiix), X e (0, L) . (4.10) 

Comparing (|4.3l)(|4.8p and (|4.10p . the initial condition in (|2.6|) implies Afc(O) = 0. Moreover, 
injecting (|4.3p(|4.8|) in the transient Schrodinger equation (|2.6p and using equations (|4.4p . (|4.5p and 
(|4.9p we get the following equation on Xk{t): 

[thX'^it) ~ z{t)Xk{t)]u{t,x) + thXk{t)dtu{t,x) = -i.ol[a3,&3](a;)^r(i>a:). 

Multiplying the previous equation by u(t,x) and integrating on (0,i), it follows: 

f Kit) + [Ht) + Io9tuit,x)uit,x)dx]Xkit)^ Skit) 
I A,(0) = 0, 

where Skit) = j-vo fl^^ ■^Y{t,x)u{t,x)dx. 

Remark 4.2. Equation (|4.1ip is an ODE which homogeneous solution oscillates at the energy 
E{t) and which source term oscillates at the energy £fc(t) defined by: 

f f^, fc>0 

Therefore Xk(t) has a peak at the frequencies k such that £fc(i) = E(t). 



4.2.1 Step 5*1 

The aim here is to compute the wave function at step I + 1 making the decomposition: 

However, at this point of the algorithm, we do not have the value of the potential y'+^ and the 
resonant mode v!"^^ can not be computed. Therefore, we will make the following approximation: 

*'+! = (*r)'+^ + ^fe"'+^ + 4:^^"'+^ ■ (4.13) 

Then the step 51 writes as follows. 

Suppose that the quantities (^J!"")', w' and Aj^, are known (at Z = 0, it is given by the initial 

decomposition (|4.10p '). 

As it is done in the Direct Resolution for ^,1,^^, see [TB], the function (^JJ*^)'^^ is computed on [0, L\ 

using (^JJ'^)' and T/'+2 by the resolution of (|4.4p with the Crank-Nicolson scheme (|B.2P , where 

the potential Q is equal to Vfm + V . Since the initial data $^'" is not supported in (0, L), the 
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suitable boundary conditions are the non-homogeneous discrete transparent boundary conditions 
presented in Appendix |B] More precisely, if the bias is equal to 

( Bi, i = 
Bit)^\ ^' , (4.14) 

[ Boo, t > 

then we will use the boundary conditions (jB.8p(jB.7p and for general biases B{t) we will use the 
boundary conditions (|B.8I)(|B.9I) . 

The function w'+^ is computed by the Cranck Nicolson scheme (jB.2|) on [0, L] starting from u' and 
using for l/'+2 the solution to (j4.ip . The initial data uj is almost equal to at x = and x ~ L 



and the homogeneous boundary conditions (|B.3p(|B.4p can be applied directly. 
The resonance z'+J and resonant mode u'+2 arc computed using the method presented in section 
13.21 where the potential Q is equal to C/'+a + l/'+2 . The initial guess for (m'+2 ^ 2'+ 2) are the first 
eigenfunction and eigenenergy of the Hamiltonian 

with homogeneous Dirichlet boundary conditions at 02 and 62. 

To achieve the step SI, we have to compute the coefficient A|^^^. It is obtained from Aj^. (^^'')', 
^^nry+i^ ^/+2 g^j^jj u'+2 by the resolution of equation (j4.11[) . The resolution is performed with a 
Crank-Nicolson scheme, which leads to the following iteration: 



(4+' - Al.)/Ai + [-z'+3 + / dtu'+-2u'+Ux]iXi + X[+^)/2 
" Jo 



&3 



{s[^^' + s[^^-'^')/2, (4.15) 



where S^^"" = ^wq / {^^^^'yT^dx. 



By adequately fixing the resonant mode phase, the quantity 

/i'+i/2:= r dtu'+^^^'^I^^^dx 
Jo 

appearing in (|4.15p is fitted to zero. Indeed, we consider u{t) solution to (j4.5[) such that ||w(t)||L2(Q j;,) : 

1 and we look for u{t) of the form u{t) = {1(^)6*"^^*' where ip{t) S M. We note first that we have 

the approximation 

2At Jq 
which becomes 



/i'+i/2 = J-Im[ / u'+^'^ui-y^dx] (4.16) 

^i Jo 



under the condition ||u' "^^^1^2(0,^) ~ I1'"'^"^^^I1l2(o,l) ~ 1- Then, defining 
we choose 



p+i72| 



and it follows: 



Jo 

As a consequence, Im[L v}^^/'^u'-~^/'^dx\ ~ and equation (j4.16p shows that u'+^Z^ is such that 
^;+i/2 jg almost equal to 0. 
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4.2.2 Step 5*2 

In this section, wc use the decomposition (j4.13p to find an approximation of 



J+i 



^M 2 

-kM 



dk 



which is adapted to the resonant peaks. Like in section [3.31 we make the approximation: 



\^ 



i+i\ 



(*r) 



i+i 



ikV 



\i+'u'+i 



2Rcid,v'+^X[+\'+-2 



(4.17) 



The coefficient 6k has a peak at Ej given by formula p.7p and as noted in Remark [4.2l the coefficient 
Afc(i) has a peak at E{t). It follows that the cross term in the previous equation can be neglected 
for biases of the form (|4.14p with Bj ^ Boo since the initial resonance and the resonance at time 
t are not in the same domaine of frequency. We write here the step S2 in the corresponding 
framework (the method for general biases follows the same line). Therefore, equation (|4.17p writes 



K^'\ 



in 



\i+i 



'^kV 



l+l\ 



\'+^u'+-2 



The wave function '^^^{t) being regular with respect to the frequency fc, the non resonant part 
of the density can be computed with a large frequency mesh size, as it is done in the stationary 
regime in [5]. Suppose that the discretisation fco = — /ca/, fci, ..., fcp, ..., fcp = /cm of the interval 
[—/cm, kM] is uniform with mesh size A/c and suppose that there exists two integers P' and v such 
that -pj ~ v. Then, the approximation for the density is 




+ 9[K(p'+i)) 



{nWr 



i+i 

3 



uAk 



A/c 



E (5(Ml^pl'+.9(fcp+i)I^P+i|') — l^ri 



p=0 

p-i 



E(5(MiA;,+^r+3(fcp+i)iAtVir)^ 



p=0 



l+h 



The number of Schrodinger equations to solve is reduced: P' equations instead of P equations for 
the Direct Resolution, which induces a reduction of the numerical cost. However, this reduction 
implies that we have only access to the functions (^"p/) for p' = 0, ..., P' and the computation 
of the coefficient A'"^^, for p = 0, ..., P, requires an interpolation of the non resonant wave function 
to evaluate the source term in (|4.15p . Following Remark 14.21 the picked form of the coefficent A'+^ 
is obtained numericaly only if the approximation of the source term at the frequency k osillates in 
time at the energy Sk(t). This implies that a polynomial interpolation is not adapted. Therefore, 
we propose the suitable algorithm below: from the scheme (|4.15p . we have for p ~ 1, ..., P 



W+i 



1 + i 



2h 



1-i- 



Atz'+i 



2h 



At + ^(^-' 



gi+^,i+i. 



where the source term is given by the interpolation: 



for < y < P' - 1 and 1 < j < i^ - I, 5^'™+^- = -vo / (^:;;0"e"^^"p'+^*"u'da;, 



63 



with 



nr\m,^j-e?°t'' 



(^„r-)m ^ (^»'-)™eTrS 
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and 




fcp > 



fco < 0. 



Remark 4.3. It is numerically verified that the method presented in this section is stable in time 
and allows long time simulations, see section [531 It is not the case of the time dependent version 
of the method presented in 13.31 which corresponds to the explicitely computation of the integral 



Ik g{k)\Xkit)\'^dk with the Briet-Wigncr formula for Afc(i) of the form (jXTt - 



5 Results 



The physical parameters used for the numerical computations are gathered in the following array: 



Rel. el. mass 


0.067 


Rel. permitivity 


11.44 


Temperature 


300 if 


Fermi level Ep 


6, 7097 X 10-21 J 


Donor density n]j 


1024 m-3 






Donor density njy 


5 X 10^1 m-3 







For all the tests, the two barriers has the same size which is equal to the size of the well. The data 
concerning the external potential arc gathered in the following array: 



L (nm) oi (nm) 02 (nm) 03 (nm) 63 [nm 
135 50 60 65 70 



,) 62 (nm) bi (nm) vq (eV) 
85 (TS 



75 



In all the simulations, we took the number of space points to be J = 300 which is such that the 
^^^ > I is verified. And we fixed ku = J^iEp + TkBT). 



stability condition 



5.1 Computation of resonances 

In this section we give the numerical values of the resonance of lower energy obtained using the 
algorithm presented in section lHT^ for different biases Bj. The potential Q is equal to Uj + Vi where 
Ui is the external potential and Vi the solution to the Schrodinger-Poisson system (|2.ip (|2.2[) (j2.4[) 
corresponding to _B/. The number of iterations before convergence Ncv denotes the iteration n 
such that ||A'/(z")u"||2 < 10~^^. We obtain the following results: 



BiieV) 


iVc. 


So(meV) 


£;7(meV) 


Ti/Ei 





3 


126.83 


127.55 


2.58 X 10-^ 


0.1 


3 


80.29 


81.00 


4.40 X 10"^ 



where the resonance is equal to zj 
Dirichlet Hamiltonian p. 8 
is represented in Figure [31 



E 



I — i^ and Eq denotes the fondamental energy of the 



The modulus of the normalized resonant mode v!" V for Br 



0.1 eF 



5.2 The stationary regime 

We show here a comparison of the Direct Resolution and of our One Mode Approximation algorithm 
presented in section [31 with respect to a Reference Resolution. The latest corresponds to the Direct 
Resolution with a uniformly refined frequency mesh where P = 4000. 
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Resonant Mode for B,=0.1 eV 



0.35 
0.3 

0.25 
0.2 

0.15 
0.1 

0.05 



-0.05 

-0.1 



N 

I ^ 

1 ' N 




:^ 


U,+V| (eV) 


1 1 ' ^ 

1 A 


- 




\ 

\ 
1 1 1 1 "~" — i_ 1 





0.2 0.4 0.6 0.8 1 1.2 

x(m) 



1.4 



xlO" 



Figure 3: Representation of the potential Ui + Vj (dashed hne) and the corresponding normalized 
resonant mode (full line) for Bj = 0.1 eV and J = 300. 



Remark 5.1. For the bias B/ = the Gummel method can be initialized at the potential V^ ~ 0. 
Such an initialization does not converge for Bi > 0, we have convergence when initializing at the 
solution given by the method when Bj = 0. 

For the One Mode Approximation, the method to compute the resonance at the first iteration 
of the Gummel algorithm is initialized at the fondamental energy and fondamental mode of the 
Dirichlet Hamiltonian (|3.8p . For the following iteration, it can be initialized at the resonance and 
resonant mode of the previous iteration to decrease the numerical cost. 

The results we obtained for two different values of the bias Bj are given in the array of Figure 
m The number of iterations before convergence Ncv denotes the iteration I such that 



W^ -V 



i-n 



\\V% 



< 10" 



(5.1) 



The integer P is the number of frequency points, CPU is the time spent by the processor to realize 



the computation, 
given by: 



The column L2 error denotes the relative error for the potential in L norm 



100- 



w 



Na. 



K-e/l 



Vr 



ref\\2 



where Ke/ is the potential of the Reference Resolution. In the Direct Resolution, the number of 
frequency points changes from an iteration to the other, however it stays around a fixed number 
which we wrote with the symbol w. The One Mode Approximation algorithm converges and needs 
less frequency points than the Direct Resolution. Nevertheless, the first algorithm doesn't cost 
much less than the second because the computation cost of one resonance is much bigger than 
the resolution of one Schrodinger equation. In the time-dependent case, it is not possible to use 
an adaptative frequency mesh and the numerical cost for the computation of one resonance is 
comparable to the resolution of one Schrodinger equation, therefore the One Mode Approximation 
will be much more interesting. 
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Ncv 


P 


CPU(s) 


L2 error(%) 


Bi 


= OeV 


Reference 


36 


4000 


69.47 


/ 


Direct 


37 


«916 


16.47 


4.59 X lO-'^ 


One Mode 


37 


50 


5.53 


1.26 


Bi = 


= 0.1 eV 


Reference 


34 


4000 


65.10 


/ 


Direct 


34 


w896 


14.80 


3.03 X 10-^ 


One Mode 


34 


50 


4.34 


2.04 



Figure 4: Comparison of the Direct Resolution and One Mode Approximation algorithms for the 
resolution of the stationary Schrodinger-Poisson system. 



40 



35 



30 



25 



20 - 



Convergence of the Gummel method 



I 15 



10 



5^ 



10"' 











One-Mode Approximation 




; 



10 



10-'" 
Error 



10 



10 



Figure 5: Relative error with respect to the number of iterations for Bi = QeV. 



In Figure [S] is represented the error e' defined by (j5.ip with respect to the number of iterations / 
for the Direct Resolution and the One Mode Approximation. The bias is equal to 5/ = eV which 
corresponds to the upper part of the array of Figure |4l We see that the One Mode Approximation 
algorithm reaches the accuracy 10~^^ in 37 iterations and converges as fast as the Direct Resolution. 
In Figure [HI are represented the self-consistant potential and the density after convergence for the 
part Bj = eV of the array of Figure |3] and for the three methods. We see that the three methods 
give very similar results. 



5.3 The transient regime: step 5*1 

We present here a validation of the step SI discribed in section lT^ It is realized using a comparison 
of the One Mode Approximation algorithm described in this section, where only the step SI is 
performed, with the Direct Resolution. For this One Mode Approximation algorithm, the integrals 
P ■ are computed with a trapezoidal rule instead of using the step S2 and the number of frequency 
points has to be the same than for the Direct Resolution. 
The external potential U{t) is the time dependent external potential corresponding to the bias 
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Comparison of the nonlinear potentials 




Companson ot the densities 




Figure 6: Nonlinear potential and density for Bj ~ OeV. 



B{t) defined in (|il^ where Bj =OeV and S^o = 0.1 eV. 

For the two methods, the number of frequency points is equal to P 



1500. In order to have 



stability, we took a time step At — 10~^^ s and the final time of the simulation is t = 8 x 10~^^ s. 
The time spent by the processor to complete the simulation is 66853.45 s for the Direct Resolution 
and 97329.36 s for the One Mode Approximation. The decomposition of the wave function at 
each frequency and time step makes the computational time more important for the One Mode 
Approximation, therefore, this method is not interesting without a reduction of the number P. 
The left plot of Figure [7] is the representation, for the two methods, of the evolution of the charge 
in the well: 



n{t, x)dx 



(5.2) 



with respect to time t. We see that for the One Mode Approximation and for the Direct Resolution 
the charge in the well increases from the charge corresponding to Bj to the one corresponding to 

Boa- 

On the right plot of Figure[71 is represented the relative distance in L2 norm d} between the density 



and the density rirx, corresponding to the stationary solution with the bias Bq 

,||'^'-"oo||L2(a2,fc2) 



100- 



jlk^laa 



(5.3) 



with respect to time t' . The relevant variation of the density being in the well, the norm is con- 
sidered only on the interval (a2, 62). We remark that for both methods, the density converges to a 
density close to the stationary density corresponding to Boo- The difference between the asymp- 
totic density and the stationary density corresponding to B^o is due to the fact that the methods 
used to perform the integration in k for the time dependent and for the stationary algorithms are 
different. 
In the left plot of Figure |8l are represented, at time t = 3 x 10^^^ s, the functions 



Ce{t,k) 




'kv{t,x)\^dx 



C\{t,k) 




Xk{t)u{t, x)\'^dx 



(5.4) 



the function C{t, k) defined by (|4.6p and the frequencies k~^[t), k'^[t) defined by (|4.7p with respect 

to the frequency k. This picture confirms that the step 51 corresponds to a decomposition of the 

wave function $/; in a function OkV peaked at the frequency related to the initial resonant energy 

El and a function XkU peaked at the frequency related to the resonant energy E[t). 

In the right plot of Figure El is represented the square of the L2 norm of v{t) in the well, defined 

by: 

N{t)^ [ ^ \vit,x)\^dx, 
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Evolution of the charge in tho weli 



Evoiulion of tiie L differenoe 




12 3 4 5 6 7 

t(s) 




Figure 7: Evolution of the charge in the well (left) and of the L^ difference d} (right) for the Direct 
Resolution and One Mode Approximation algorithms. 



Ctiarge in ffie weii wiffi respeot to ffie frequenoy af fime t=3x 10 ^^s 




— C(f,li) 

- -G,(t,k) 

Cj(t,k) 
- *r('> 




Figure 8: Functions C(i, fc), Cg{t,k), C\{t,k) and resonant frequencies kj^{t), k^{t) at time t 
3 X 10"^^ s (left) and time evolution of N{t) (right). 



as a function of the time t. We observe that v{t) vanishes in the well for times larger than 
5 X 10~^^ s. This corresponds to the extinction of the resonant peak at energy Ej observed in 
Figure [21 Moreover, the decay of N{t) is theoritically given by the imaginary part of the resonance 
following the approximation: 

N{to) 
This is verified numerically since we have for to = 10~^^ s and T = 10~^^ s: 



Njto + T) 
Nito) 



0.578 



and 



e-*^r"r(.)<:^^ = 0.582. 



5.4 The transient regime: reduction of the number of frequency points 



We present here a comparison of the One Mode Approximation algorithm described in section 14^ 
where the steps SI and S2 are performed, with the Direct Resolution. 
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Evolution of Ihe charge in thG well 




Evolution of the L differencG 




Figure 9: Evolution of the charge in the weh (left) and of the L^ difference d} (right) for the Direct 
Resolution and One Mode Approximation algorithms. 



The external potential U{t) is the time dependent external potential corresponding to the bias B{t) 
defined in (|4.14[) where Bj = eV, Boo =0.1 eV. The number of frequency points is P = 1500 for 
the Direct Resolution. For the One Mode Approximation, the number of frequency points for the 
non resonant part is P' = 750 and P — 1500 for the resonant part. As in the previous section, we 
took a time step At ~ 10~^^ s to have stability. The final time of the simulation is t = 8 x 10~^^ s. 
The left plot of Figure IHl is the representation, for both methods, of the charge evolution in the 
well defined by (|5.2p . We see that for the One Mode Approximation and for the Direct Resolution 
the charge in the well increases from the charge corresponding to Bj to this corresponding to Boo ■ 
On the right plot of Figure [HI is represented the L2 difference d', defined by (|5.3p . as a function 
of time tK We remark that for both methods, the density converges to the stationary density 
corresponding to Boo ■ 

For the One Mode Approximation the simulation is possible with a half less Schrodinger equations 
at each time step and the computational time is divided by two. The method we propose for the 
step 5*2 presents the two following properties: first, the resolution is stable and allows long time 
simulations, second, the picked form of the coefficient is reconstructed and allows at time increase 
of the charge in the well. These two properties, that we set as a validity criterion for the One Mode 
Approximation algorithm, are still verified after a reduction of the frequency points similar to the 
stationary algorithm by taking P' = 50 and P = 1500. However, for such a simulation, the intial 
value of the density and the long time value of the density are poorly evaluated. Nevertheless, this 
lack of precision can be probably overcomed by improving the interpolation of the exterior wave 
functions and replacing the trapezoidal rule to compute the initial resonant pic J^^'''^^ g{k)\9k\'^dk 
by its explicit value given in section 13.31 

We represented in Figure [TU] the logarithm C{t, k) of the charge at frequency k, the frequencies 
k~^{t) and k^{t) defined by (j4.6p and (j4.7p computed with the Direct Resolution and at different 
values of the time t. As noted for Figure [21 for small times the peak is localized at the initial 
resonant energy Ej. Then, the initial resonance peak vanishes while a resonance peak grows 
around E{t). However, the two cases are different: in Figure [5] the resonant energy E{t) moves 
with the bias B{t) for t £ (0,to)- In that case, the resonance peak corresponding to E{t) appears 
at the same frequency as the initial one and the two peaks split. In the present case, the peak, 
when it appears, is already far from the initial one. 
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Figure 10: Logarithm C(t, fc) of the charge in the weh for one wave function with respect to the 
frequency fc and frequencies, kjj{t) and k]^{t), corresponding to the resonant energy at different 
time t. 
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A Computation of resonances: FEM matrices 

The matrices appearing in section [3.21 are given as follows: 

\ 



A/2 = 



-t^^ ... \ f 

M3 = 



I ~ \/2m 





V 







) 











Ma = Ax 



( 1 -1 

-1 2 -1 



Ml 



IvtiAx 



V ... -i^ / 



/ Co Co 

Co ix Ci 



/ ? i 

6 3 





6 5 5 

6 3 / 



-1 2 -1 

-1 1 / 







Cj-2 O-l C./-1 

Cj-i O / 



where 



Qj + Qj+i ,• - J - 1 

12 ' ^-"'-'"^ ^' 



Qo , Qi ^ Qj-i + Qj+i , Qj , ^ 7 1 ^ 0./ , Q./-1 



^''"T + T2' ^' 
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-,, = 1,...,J-1, = —, ^2 



and Qj are the approximations of the nodal values of the potential: Q{xj). 

B Resolution of the time dependent Schrodinger equation 

B.l The homogeneous case 

In this section, we recall the scheme proposed in |2 to solve on the bounded domain [0,L] a time 
dependent Schrodinger equation 



(B.l) 



( ihdt^ = -^dl^ + Q^, t>0,xeM 
\ *(0,a;) = $(a;), x e M, 

with the following hypothesis: 

HI. The initial condition $ is supported in < x < L. 
H2. The potential Q verifies: for f > 

Q{t, x) — 0, X <0, and Q{t, x) ~ Ql, x > L . 

Considering the time and space discretization defined in section HTTl we note ^' the approximation 
of the solution 'i>{t'',Xj). Then, equation (jB.l|) . is solved with the Crank- Nicolson method: 



' + l/2/,T,' + l_L,T,'- 



iR{¥+' - *;.) - A,,*^+^ + A,,*^. + u;Q7^/^(*^+^ + *;.) 



j = l,...,J-l, />0, (B.2) 



where A,,*^ = *j+i - 2*, + *,_i, R = ^^^ and w = -^^^. 

The equation (|B.2[) comes with the discrete transparent boundary conditions: 

(-1 



*'i-5o< = E4''*o-*'r\ ^>i, 



(B.3) 



k=l 



21 



l-l 



* 



,/-l 



'"*S = E-^'r'* 



k^jfk 



* 



j-i' 



/>i, 



where, we have for j = and j ^ L 

R a 



1 -I h — 



<5P + 



1 + i— + — 
2 2 



fe=i 



^.^ 



Qfj CXp 



-Hip 



Flip,) - Pl-2i^i,) 



21- 1 



and 



ipj = arctan 



2R{aj + 2) 



i?2 + 4cr,- + cr? 



Mj 



2Aa 



h^ 



-Q, 



^iR^+a])[R^ + ia,+Ar] 
{{R' + a^)[R' + {a, + if])'/' o.p {^^) 



(B.4) 



(B.5) 



0, and S-j the 



Here Pi denotes the Legcndre polynomials, with the convention P_i = P_2 

Kronecker symbol related to the integers j and I. 

The boundary conditions above are available only for initial data supported in (0, L) and therefore 

it is called homogeneous discrete transparent boundary conditions. 



B.2 The non homogeneous case 

For the intial potential Qj such that Q(0, x) = Qi{x) and 

Qi{x) = 0, a; < 0, and Qi{x) = Qi.l, x> L , 
we consider the problem (jB.ip where the initial condition $ is solution to 

-^— -— $ + Q/$ = S7$, xelR., 
Zm dx'^ 



(B.6) 



for a given energy i?/. Like in section IBTTI we will use the Crank-Nicolson scheme (|B.2|) . however, 
the hypothesis ffl is not verified and we can not apply (|B.3p(|B.4[) . 
Nevertheless, using equations (jB.ip and ()B.6|) . the function 



is solution to: 



V? = * - $e '— * where El = Ej + {Ql - Qi,l) 

ihdttp = h7;—dl + Q]ip, x>L 
Zm 



and verifies (p(0, x) = 0. Thus, we can write the homogeneous boundary condition (jB.41) for ip 
which gives the following boundary condition at a; = L for ^: 



i-i 



fe=l k=l 

+ $j_ie-*^('-i''^*(l + e-'^^*), 1>1, (B.7) 
where the coefficients s' are given by (|B.5p . We proceed similarly at x = by setting 

(^ = * - $£ ^ h ^ where Eq ^ Ej . 
We obtain the following boundary condition at a: = 0: 



fc=i 



fc=i 



+ <I>ie-'T^('-i)^*(l + e-*^^*), />1. (B 
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B.3 Time dependent potential 

In this section, we consider the same problem than in section [B]2] and in addition, we suppose that 
the potential verifies: for i > 

Q{t,x) =0,x<0, and Q{t,x) ^ QL{t), x> L , 



instead of the hypothesis H2. As we did in section |B.2[ we will use the Crank-Nicolson scheme 
(|B.2[) . The boundary condition (jB.8|) at x = obtained there can also be applied in our case. At 
X = L, the time depcndancc of the potential must be considered. To this aim, we introduce the 
function: 

(^ = ^e*^o'5^("'''^-$e-'-^* where El^Ej-Qi.l. 

Using equations (jB.ip and (jB.6|) . it is verified that (^(0) =0 and 

ihdtif = -T^—dlif, x> L , 
Zm 

where we got rid of the time dependent potential. Therefore, ip verifies at a; = L the boundary 
condition (|B.4p with potential equal to 0. This gives the following boundary condition for ^: 

yi > 1 



/-I 



Sj e B 



fe=l k=l 

+ $j_ie-'^('-i)^'(l + e-'*^*) , (B.9) 
where 

and s''j is given by (jB.51) when replacing ctj by 0. 

C Explicit value of some useful integrals 

The aim of this section is to give, for n = 0, 1, the explicit value of the integral 

f'' x" 

where a < b, c > and d > 0. 



• Case n=0: 

We have the decomposition 



1 -1/1 1 



(x^ — c)^ + d^ 2id \x'^ — c + id x"^ — c — id 



^Im ^ ' 



and 



d Kx"^ — c + id 



1 -1/1 1 



x'^ ~ c + id 2zq \x + zq X — zq 
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(C.l) 



where zq — Vc — id and ^/z denotes the square root holomorphic on C \ IR- . It is defined 
by z = pe'^ and -y/i = \fpe}^ where p > 0, € (— tt, tt]. 
This leads to: 



X°(a,6,c,d) = — Ini 



2d 



I f fb I fb I 



dx — dx 



^0 \Ja X + Zq J ^ X — Zq 

— (ln(& + zo) - ln(a + zq) - ln(6 — zq) + hi(a — zq)) 
^0 



where the logarithm is defined by z = pe'^ and \x\z = \n p + iO where p > 0, G (^""j ^]- 

• Case n=l: 
We have 

,2 

1 /■ 1 

2 Ja2-c x^ +d^ 

1 / /&'-c\ /a2_^ 

-— arctan | ; — | — arctan 



2d V \ d 
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